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Carnegie  Mellon  University  is  automating  the  use  of  Ground  Penetrating  Radar  (GPR)  for  the  cleanup  of  hazardous  waste 
sites.  Tlic  current  focus  is  on  the  dcvelopmem  of  an  automated  subsurface  mapping  system  to  locate  buried  objects  and 
geological  structures  so  that  sources  and  migratory  pathways  of  contaminants  can  be  identified  and  cataloged.  The 
subsurface  maps  are  produeed  using  the  non-invasive  sensing  capabilities  of  GPR.  GPR  operates  similar  to  conventional 
radar,  but  the  data  acquired  is  more  difficult  to  process  due  to  the  nonhomogencous  nature  of  the  subsurface  medium. 
Though  sometimes  used  in  waste  site  characterization,  GPR  deployment,  data  aequisition,  and  interpretation  are  human 
driven  processes.  The  potential  of  GPR  to  generate  accurate  three  dimensional  subsurface  maps  has  not  been  fully 
realized  previously.  The  Site  Investigation  Robot  uses  robots  to  position  a  GPR  transducer  to  exploit  the  accurate, 
repeatable  positioning  available  from  automated  equipment.  By  combining  the  use  of  a  position  cognizant,  all  terrain 
mobile  robot  and  a  linear  scanning  mechanism,  GPR  records  arc  acquired  in  a  two  dimensional  grid  on  the  ground 
surface.  This  method  of  collection  simplifies  processing  and  the  positional  location  of  features  of  interest  in  the  data. 

To  achieve  accurate  positional  accuracy  of  the  located  subsurface  objects,  correct  modelling  of  die  subsurface  medium  and 
the  antenna  is  crucial.  The  subsurface  medium  can  be  modeled  as  a  composition  of  multiple  subvolumes,  each  one  being 
characterized  by  its  own  electrical  parameters.  The  radar  wave's  energy  can  be  calculated  at  any  ptTint  in  the  subsurface 
using  these  parameters  and  tlic  distance  the  wave  traveled.  Analysis  of  die  wave's  propagation  in  die  subsurface  is 
simplified  by  the  u.se  of  a  tcchnii|uc  known  as  ray  tracing.  Two  criteria  arc  used  to  determine  if  die  u.sc  of  ray  tracing  is 
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valid:  the  wavelength  of  the  radar  waves  must  be  short  in  comparison  to  the  resolution  of  the  external  perturbations  of 
any  expected  objects  and  the  medium  must  be  an  isolator.  Modelling  of  the  GPR  antenna  is  best  inmxJuccd  by  the 
analogy  of  a  laser  and  a  flash  light. 
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Abstract 


Carnegie  Mellon  University  is  automating  the  use  of  Ground  Penetrating  Radar  (GPR)  for  the  cleanup  of  hazardous 
waste  sites.  The  current  focus  is  on  the  development  of  an  automated  subsurface  mapping  system  to  locate  buried 
objects  and  geological  structures  so  that  sources  and  migratory  pathways  of  contaminants  can  be  identified  and  cata¬ 
logued.  The  subsurface  maps  are  produced  using  the  non-invasive  sensing  capabilities  of  GPR.  GPR  operates  similar 
to  conventional  radar,  but  the  data  acquired  is  more  difficult  to  process  due  to  the  nonhomogeneous  nature  of  the  sub¬ 
surface  medium.  Though  sometimes  used  in  waste  site  characterization,  GPR  deployment,  data  acquisition,  and  inter¬ 
pretation  are  human  driven  processes.  The  potential  of  GPR  to  generate  accurate  three  dimensional  subsurface  maps 
has  not  been  fully  realized  previously.  The  Site  Investigation  Robot  uses  robots  to  position  a  GPR  transducer  to 
exploit  the  accurate,  repeatable  positioning  available  firom  automated  equipment.  By  combining  the  use  of  a  position 
cognizant,  all  terrain  mobile  robot  and  a  linear  scanning  mechanism,  GPR  records  are  acquired  in  a  two  dimensional 
grid  on  the  ground  surface.  This  method  of  collection  simplifies  processing  and  the  positional  location  of  features  of 
interest  in  the  data.  To  achieve  accurate  positional  accuracy  of  the  located  subsurface  objects,  correct  modelling  of  the 
subsurface  medium  and  the  antenna  is  crucial.  The  subsurface  medium  can  be  modelled  as  a  composition  of  multiple 
subvolumes,  each  one  being  characterized  by  its  own  electrical  parameters.  JQie  radar  wave’s  energy  can  be  calcu¬ 
lated  at  any  point  in  the  subsurface  using  these  parameters  and  the  distance  thdi  wave  travelled.  Analysis  of  the  wave’s 
propagation  in  tlte  subsurface  is  simplified  by  the  use  of  a  technique  known  as  ray  tracing.  Two  criteria  are  used  to 
determine  if  the  use  of  ray  tracing  is  valid:  the  wavelength  of  the  radar  waves  must  be  short  in  comparison  to  the  res¬ 
olution  of  the  external  perturbations  of  any  expected  objects  and  the  medium  must  be  an  isolator.  Modelling  of  the 
GPR  antenna  is  best  introduced  by  the  analogy  of  a  laser  and  a  flash  lighL  A  laser  is  a  highly  focused  beam  of  light, 
but  tlie  intensity  of  the  light  transmitted  by  a  flash  light  decreases  rapidly  with  distance  from  the  center  of  the  light 
The  GPR  antenna  (grates  in  a  similar  fashion  to  the  flash  light  Two  models  of  the  antenna  are  used:  the  energy  of 
the  antenna’s  beam  is  assumed  constant  at  all  points  equidistant  from  the  antenna  and  the  empirical  measurements  of 
the  beam’s  attenuation  pattern.  Because  of  the  lack  of  focus  in  the  antenna’s  beam,  the  data  collected  at  one  surface 
location  represents  reflections  from  all  objects  and  interfaces  within  the  beam.  Migration  is  the  signal  processing 
operation  used  to  reverse  a  spatial  integration  effect  in  the  collected  data.  Homogeneous  migration  assumes  the  elec¬ 
trical  characteristics  of  the  subsurface  are  constant  throughout  the  entire  volume  processed;  heterogeneous  migration, 
a  more  realistic  operation,  allows  these  parameters  to  vary.  Toward  Modelling  simulates,  on  a  computer,  the  collec¬ 
tion  of  radar  data  under  user  specified  conditions.  There  are  two  purposes  for  forward  modelling:  test  the  migration 
operation  and  as  a  learning  tool  for  understanding  the  propagation  of  GPR  waves  in  the  subsurface.  Because  of  the 
design  of  the  software,  the  user  can  easily  change  parameters  that  control  the  simulation  (e.g.  -  antenna  height).  The 
accurate  tracing  of  radar  waves  requites  a  model  which  includes  the  reflection  and  refraction  of  rays  at  subsurface 
interfaces  and  objects.  This  model  is  a  necessity  for  heterogeneous  migration  and  forward  modelling  (not  currently 
implemented). 


1 .0  Executive  Summary 


1.1  statement  of  Purpose 

Carnegie  Mellon  University  is  automating  the  use  of  Ground  Penetrating  Radar  (GPR)  for  the  cleanup  of  hazardous 
waste  sites.  The  current  focus  is  on  the  development  of  an  automated  subsurface  mapping  system  to  locate  buried 
objects  and  geological  structures  so  that  sources  and  migratory  pathways  of  contaminants  can  be  identified  and  cata¬ 
loged.  The  subsurface  maps  produced  use  the  non-invasive  sensing  capabilities  of  GPR.  GPR  operates  on  principles 
similar  to  conventional  radar,  but  the  data  acquired  is  more  difficult  u>  process  due  to  the  heterogeneous  nauire  of  the 
subsurface  medium.  Though  sometimes  used  in  waste  site  characterization,  GPR  deployment,  data  acquisition  and 
interpretation  are  human  driven  processe.s.  The  potential  of  GPR  to  generate  accurate  three  dimensional  subsurface 
maps  is  thus  not  fully  realized  in  practice.  The  Site  Investigation  Robot  project  uses  robots  to  position  GPR  transduc¬ 
ers  to  exploit  the  accurate,  repeatable  positioning  available  from  automated  equipment.  By  combining  the  use  of  a 
position  cognizant,  all-terrain  mobile  robot  and  a  linear  scanning  mechanism,  it  is  possible  ro  acquire  GPR  records  in 
a  two  dimensional  grid  on  the  ground  surface. 

Data  collected  using  the  GPR  anteima  does  not  accurately  reflect  the  structure  of  the  subsurface  for  two  reasons.  The 
volumetric  nature  of  the  antenna  beam  causes  spatial  integration  of  the  subsurface  structure  before  the  data  is  col¬ 
lected.  The  second  reason  for  inaccuracy  in  the  reflected  data  is  distortion.  During  propagation,  the  radar  wave  is  con¬ 
volved  with  the  linear  system  which  represents  the  structure  of  the  subsurface.  Also,  the  antenna  pulse  is  not  a  perfect 
impulse  function;  this  causes  the  collected  data  to  be  convolved  with  imperfections  in  the  antenna  pulse.  Both  of 
these  problems  are  addressed  by  the  deconvolution  operation  (not  explained  in  this  paper). 

Two  GPR  antennas  are  used:  one  for  transmission  and  one  for  reception  of  radar  waves.  High  frequency  radar  waves 
transmitted  by  the  send  antenna  are  reflected  back  to  the  receive  antenna  by  interfaces  between  subsurface  materials 
with  different  electromagnetic  characteristics.  Reflections  received  from  the  same  point  on  the  surface  are  stored  in  a 
one  dimensional  data  array,  and  tagged  with  the  surface  position.  All  the  GPR  data  collected  over  the  entire  two 
dimensional  grid  form  a  volume.  The  operation  of  the  GPR  antenna  can  be  illustrated  using  a  laser  and  a  flashlighL  A 
laser  is  a  highly  focused  beam  of  light,  but  the  light  intensity  of  a  flashlight  decreases  rapidly  with  distance  from  the 
center  of  the  beam.  The  GPR  antenna  operates  in  an  similar  fashion  to  the  flashlight.  Because  of  this  lack  of  focus  in 
the  antenna  beam,  the  data  collected  at  one  grid  location  represents  reflections  from  all  objects  and  interfaces  within 
the  antenna’s  beam.  Migration  is  the  signal  processing  operation  used  to  reverse  this  distortion  in  the  collected  GPR 
data.  Due  to  the  heterogeneous  nature  of  the  subsurface,  the  exact  shape  of  buried  objects  can  never  be  fully  recov¬ 
ered  by  the  migration  operation.  Conventional  image  processing  operations  are  used  on  the  post  migrated  data  to 
decrease  the  needed  for  a  highly  skilled  interpreter.  Forward  Modelling  simulates,  on  a  computer,  the  collection  of 
radar  data  under  user  specified  conditions.  There  are  two  purposes  for  modelling  GPR:  test  the  migration  operation 
and  aid  in  the  understanding  of  how  GPR  waves  propagate  in  the  subsurface. 

1 .2  Summary  of  Results 

Since  the  GPR  antenna’s  beam  is  three  dimensional,  three  dimensional  migration  is  the  most  appropriate  technique  to 
reverse  the  spatial  integration  phenomena.  Testing  of  the  algorithm  was  done  using  real  data  from  a  test  pit  (Figure  1) 
filled  with  dry  sand  and  a  metallic  drum  and  data  generated  from  forward  modelling. 

The  migration  operation  causes  positive  reinforcement  in  the  area  where  reflections  due  to  subsurface  objects  or  inter¬ 
faces  occurred.  However,  the  migration  operation  causes  negative  reinforcement  in  the  areas  where  there  are  no 
reflections  due  to  subsurface  objects  or  interfaces.  For  areas  where  there  are  only  weak  reflections  the  interference 
caused  by  the  migration  operation  will  produce  noise.  Real  GPR  data  (Figures  2,  3)  contains  values  with  positive  and 
negative  signs,  thus  allowing  migration’s  constructive  and  destructive  interference  to  work.  Because  the  Forward 
Modelling  was  done  in  the  time  domain,  the  generated  radar  data  only  contains  positive  signs.  Therefore,  only  the 
constructive  interference  works  when  migration  is  applied  to  data  generated  by  forward  modelling  (Figiues  4,5). 
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HGURE  1 


Configuration  of  Test  Pit 


RGURE  2 


Unprocessed  Data  front  the  Test  Pit 
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FIGURE  3 


3D  Migrated  Data  from  the  Test  Pit 
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FIGURE  4 


Unprocessed  Data  Generated  by  Forward  Modelling 
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FIGURE  5 
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1.3  Summary  of  Recommendations 


Impicmcntaiion  of  three  dimensional  heterogeneous  migration  and  forward  modelling  is  the  next  important  step. 
Since  the  velocity  of  GPR  waves  varies  significantly  in  different  medium,  the  positional  accuracy  of  the  located  sub- 
-surface  objects  should  significantly  improve  using  a  heterogeneous  algorithm.  Current  evaluation  of  the  algorithm’s 
accuracy  is  done  in  a  semi-quantitative  manner  another  important  conuibution  is  a  system  to  perform  automatic 
quantitative  evaluation  of  the  algorithm's  accuracy. 


2.0  Modelling  of  the  Subsurface  Medium 

Analy.sis  of  GPR  wave  propagation  is  done  using  two  tools.  Electromagnetic  wave  theory  is  the  tool  which  best  mod¬ 
els  the  physics  involved.  The  mathematics  are  simpler  if  ray  tracing,  a  technique  from  geometric  optics,  can  be  used. 
(Nayar  1989)  Two  cntcria  (Ulriksen  1982)  arc  used  to  justify  the  use  of  ray  uacing;  the  wavelength  of  the  radar 
waves  is  much  less  than  the  resolution  of  the  cxtcmal  perturbations  for  any  buned  object  and  the  medium  is  an  is<ila- 
tor.  Since  the  smallest  expected  wavelength  is  6  centimeters,  thc.se  two  critena  arc  satisfied. 
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FIGURE  6 


Illustration  of  Subsurface  GPR  Wave  Propagation 


GPR  Antenna 


Ray  Tracing  Justification  Criteria 


X 


c 

(/xje) 


where  f  =  500  Mhz  and  c  -  3x10*  meters  per  second 


Examples  Dry  Sand  (e  =  4  )  X  =  30  cm  Saturated  Sand  (e  =  81)  X  =  6  cm 


It  is  desired  to  trace  the  propagation  of  radar  waves  in  the  subsurface.  Since  it  is  valid  to  use  ray  tracing,  tracing  the 
propagation  of  waves  reduces  to  the  tracing  of  vectors.  Modifications  to  pure  ray  tracing  are  based  on  the  antenna’s 
beam  pattern  (see  section  3).  Figure  7  illustrates  the  principles  of  ray  uacing  (Bom  1980)  at  a  subsurface  interface. 
The  angles  of  reflection  and  refraction  are  a  function  of  the  medium’s  electrical  characteristics.  The  refracted  ray  lies 
in  the  same  plane,  called  the  plane  of  incidence,  as  the  incident  ray  and  the  normal  to  the  surface.  The  reflected  ray 
lies  in  the  plane  of  incidence. 
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RGURE  7 


Illustration  of  Ray  Tracing 
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3.0  Antenna  Modelling 

There  are  three  factors  that  affect  the  amplitude  of  the  transmitted  GPR  signal;  the  bandwidth,  the  beam  pauem,  and 
the  radar  beam’s  attenuation  in  the  subsurface  medium.  The  bandwidth  describes  tliC  change  in  signal  amplitude  ver¬ 
sus  frequency.  The  500  Mhz  antenna  has  a  bandwidth  of  ±2SQMh2.  This  means  that  the  half  power  point  (3  db)  is  at 
250  and  750  Mhz.  Since  all  the  processing  in  this  paper  was  done  in  the  time  domain,  the  bandwidth  of  die  antenna 
was  not  used. 

The  beam  pattern  describes  the  attenuation  of  the  signal  which  is  a  result  of  the  unfocused  nature  of  the  antenna  pulse. 
In  part  1,  all  energy  in  the  beam  pattern  is  assumed  to  be  constant  A  more  accurate  model  (part  2)  takes  into  account 
the  geometrical  spreading  losses  of  the  radar  wave. 

3.1  Constant  Energy  Beam  Pattern 

The  first  approximation  of  the  beam  pattern  assumes  the  radar  wave’s  energy  is  constant  for  all  points  equidistant 
from  the  antenna.  The  range  of  the  antenna’s  beam  pattern,  as  obtained  from  the  manufacturer  (Operations  Manual 
1987),  is  ±45“’  for  the  x-axis  cross  section  and  +30“’  for  the  y-axis  cross  section  (Figure  8).  The  analogy  of  the  arc  for 
the  three  dimensional  case  is  a  sphere  bounded  by  an  elliptical  cone.  One  x-y  section  taken  from  the  intersection  of 
the  sphere  with  the  elliptical  cone  is  an  ellipse. 


FIGURE  8 


Beam  Patterns  Assuming  Constant  Energy 


The  equation  for  the  aforementioned  cone  is  shown  below  (spherical  coordinate  system). 
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Since  the  radar  energy  flux  is  uniformly  distributed  across  the  surface  of  the  sphere,  accurate  three  dimensional  pro¬ 
cessing  requires  equal  weighting  of  each  ray  in  the  radar  beam.  In  order  to  accomplish  this,  the  sphere  is  divided  into 
trapezoidal  patches  of  equal  area.  Figure  9  shows  a  hemi-sphere  with  the  lengths  of  the  edges  labelled  for  a  sample 
trapezoid.  The  algorithm  divides  the  sphere’s  surface  into  a  user  specified  number  of  longitudinal  strips.  Each  strip  is 
divided  into  as  many  trapezoidal  patches  of  the  specified  area  as  possible.  The  algorithm’s  performance  was  mea¬ 
sured  by  calculating  the  ratio  of  the  area  covered  by  the  trapezoids  to  the  total  area  possible  in  the  sphere  (96.48%). 


FIGURE  9 
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3.2  Beam  Attenuation  Pattern 

Four  sets  of  beam  pattern  measurements  were  obtained:  one  from  Geophysical  Survey  Systems  Incorporated  (Shutz 
1990)  and  three  from  Steven  Duke’s  thesis  (Duke  1990).  The  measurements  obtained  from  GSSI  characterize  all  the 
GPR  antennas  (medium  of  air)  they  manufacture.  The  data  from  Duke’s  thesis  describe  the  500  Mhz  antenna  in  three 
different  medium. 

The  beam’s  attenuation  pattern  is  characterized  by  two  orthogonal  cross  sections  (Figure  10).  The  contribution  of  the 
beam  pattern  to  the  attenuation  of  the  wave  is  a  function  only  of  the  angle  from  the  z  axis  and  is  independent  of  the 
distance  from  the  antenna.  The  effect  due  to  the  distance  of  the  wave  from  the  antenna  is  described  by  the  attenuation 
equation.  Since  the  migration  algorithm  is  three  dimensional,  a  volumetric  beam  pattern  was  needed.  This  need  was 
supplied  by  using  a  volumetric  interpolation  algorithm.  The  algorithm  uses  two  steps:  a  standard  cubic  spline  interpo¬ 
lation  algorithm  and  a  custom  elliptical  interpolation  algorithm.  The  resolution  of  the  beam  pattern  cross  sections  pro¬ 
vided  was  too  low;  therefore,  the  cubic  spline  algorithm  was  used  to  increase  the  number  of  (angle,  attenuation)  pairs. 

The  crucial  insight  used  to  develop  the  custom  algorithm  is  now  described.  A  set  of  concentric  ellipses  with  centers  at 
different  values  of  attenuation  (z  axis)  can  be  described  using  the  x  and  y  axis  cross  sections  (Figure  1 1 ).  The  z  axis  is 
divided  into  points  of  equal  resolution. 


Area  of  a  Trapezoidal  Patch 

rxd^x  sin0^_, 


X  (sin0,_,  +sin0,  )  (Ax  dQ  -  d(^^  x  (sin0^_,  -  sin0j2)  ) 
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Format  of  Beam  Pattern  Measurements 
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HGURE  11 


Step  TWO  In  Volumetric  Interpolation 
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4.C  Migration 


4.1  Motivation 

Migration  is  a  processing  operation  that  reverses  the  effect  of  the  unfocused  nature  of  the  antenna  beam  in  the  col¬ 
lected  data.  In  section  2  it  was  shown  that  it  is  valid  to  use  the  ray  tracing  technique  horn  geometric  optics  in  analyz¬ 
ing  GPR  wave  propagation.  Figure  12  is  used  to  show  the  need  for  the  migration  processing  operation  (Operations 
Manual  1987);  the  example  is  two  dimensional  for  simplicity.  The  dotted  line  rectangle  represents  the  one  dimen¬ 
sional  data  array  which  contains  GPR  data  collected  from  one  position  on  the  surface.  The  data  array  is  tagged  with 
the  surface  position,  but  all  the  data  in  the  array  is  not  ?  result  of  reflections  from  that  position.  Each  entry  in  the  data 
array  represents  all  objects  and  interfaces  equidistant  from  the  antenna  Two  rays  will  be  traced  in  this  example.  Ray 
1  strikes  object  A,  and  is  reflected  back  to  the  antenna  with  a  time  delay  of  r, .  Ray  2  strikes  object  B,  and  is  reflected 
back  to  the  antenna  with  a  time  delay  of  <2-  For  homogeneous  medium,  the  distance  the  ray  travelled  is  calculated  by 
multiplying  the  time  delay  by  the  velocity  of  light  in  the  subsurface  medium  (Operations  Manual  1987).  This  distance 
is  convmed  to  an  index  into  the  data  array.  The  top  data  array  entry  contains  data  representing  objects  and  interfaces 
50  picoseconds  away  from  the  antenna.  Every  other  data  array  entry  contains  data  collected  at  an  integer  multiple  of 
SO  picoseconds  away  from  the  antenna.  As  Figure  12  shows,  reflections  of  GPR  waves  that  are  stored  in  this  data 
array  can  be  from  any  object  within  the  antenna  beam.  The  problem  just  described  is  also  true  for  the  three  dimen¬ 
sional  case. 


FIGURE  12  Motivation  for  Migration 

GPR  Antenna 


=  distance  index  corresponding  to  /] 
^2  =  distance  index  corresponding  to 


4.2  Explanation 

For  ease  of  explanation,  the  following  description  of  migration  is  done  for  the  two  dimensional  case.  Changes  needed 
for  the  three  dimensional  case  are  explained  later.  Migration  is  done  in  three  steps:  calculate  a  template,  use  the  tem¬ 
plate  to  process  the  section  of  data,  and  normalization  (Christian  1990).  One  side  of  the  array  corresponds  to  GPR 
data  collected  near  the  surface  of  the  ground.  The  template  (Figures  13)  consists  of  a  sequence  of  concentric  circular 
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arcs  of  different  radii;  the  center  is  the  data  entry  corresponding  to  the  surface  of  the  ground  and  the  radius  is  the  dis¬ 
tance  hrom  the  antenna.  Each  arc  is  subdivided  into  a  sequence  of  points;  the  distance  between  each  point  is  the  same. 
The  point  directly  below  the  data  array  element  is  called  the  source  element;  the  rest  of  the  elements  are  called  desti¬ 
nation  elements. 


RGURE  13  Explanation  of  Migration 


Side  of  Array  Corresponding  to  the 
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Review  of  Antenna  Beam  Pattern  Models 


The  first  approximation  to  the  antenna’s  beam  pattern  assumes  the  radar  wave’s  energy  is  constant  for  all  points  along 
each  arc’s  edge.  This  assumption  is  the  reason  for  the  addition  of  100%  of  each  source  element  to  all  of  its  destination 
elements.  A  more  accurate  model  of  the  antenna’s  beam  pattern  is  shown  in  Figure  13b.  In  this  model,  the  energy  of 
the  antenna’s  beam  is  at  its  maximum  value  directly  below  the  antenna.  Along  an  arc  of  constant  depth,  the  beam’s 
energy  decreases  as  the  distance  from  the  source  element  incrcases.The  attenuation  fractions  can  be  used  in  the 
migration  algorithm  for  weighting  the  source  element  before  adding  it  to  a  destination  clcment.The  extension  of 
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migration  to  three  dimensions  can  now  be  stated.  The  three  phases  of  the  algorithm  (calculate  the  template,  use  the 
template  in  processing,  nonnalization)  are  the  same.  The  template  used  is  calculated  using  one  of  the  algorithms 
described  in  section  2:  sphere  tessellation  or  volumetric  interpolation  of  the  antenna  beam’s  attenuation  pattern.  The 
attenuation  fraction  is  used  as  a  scale  factor  (saiite  method  as  two  dimensional  migration). 


5.0  Forward  Modelling 

Forward  Modelling  simulates  the  generation  of  GPR  data  on  a  computer  (Duke  1990).  Parameters  characterizing  the 
antenna,  the  subsurface,  and  the  geometry  of  buried  objects  are  easily  changed.  This  feature  makes  the  forward  mod¬ 
elling  software  a  powerful  scientific  tool.  Figure  12  was  previously  used  to  explain  the  motive  for  migration.  It  is  also 
an  example  of  how  the  three  dimensional  simulator  works.  For  flexibility  of  implementation,  the  simulation  used  a 
public  domain  ray  tracer  to  perform  the  ray-object  intersections. 


5.1  Simulation  of  Antenna  Movement 

Forward  modelling  of  GPR  data  was  done  for  a  rectangular  site.  Both  antennas  are  moved  as  a  unit  (constant  antenna 
separation  distance  -  Figure  14).  The  site  is  divided  into  a  rectangular  grid.  Simulation  oi  the  GPR  antenna’s  opera¬ 
tion  is  done  or  each  cell  in  the  grid.  The  next  section  explains  the  simulation  of  the  antenna. 


FIGURE  14  Antenna  Movement  in  Forward  Modelling 
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5.2  Reflection  Model 


All  rays  originate  at  the  receive  antenna,  which  is  modelled  as  a  point.  For  simplicity  of  implementation,  the  send 
antenna  is  modelled  as  a  circular  aperture  (Pilant  1990)  The  rays  are  traced  using  a  constant  angular  resolution  for  6 
and  <t>.  Rays  originate  at  the  receive  antenna  and  terminate  at  the  send  antenna;  consistency  with  the  original  ray 
tracer  implementation  is  obtained  by  using  this  method.  When  a  ray  intersects  an  object,  the  reflected  ray  is  calculated 
assuming  a  specular  reflection  model.  Next,  the  antenna  center  ray  is  calculated.  Its  origin  is  the  object’s  intersection 
point  and  its  end  point  is  the  center  of  the  send  antenna’s  aperture.  The  maximum  size  of  the  angle  between  the 
reflected  ray  and  the  antenna  center  ray  (3)  is  determined  by  the  size  of  the  send  antenna’s  aperture  (Shutz  1990). 


RGURE  15  Reflection  Model  for  Forward  Modelling 
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5.3  User  Created  Simulation  Scenario 

Before  making  a  simulation  run,  the  configuration  hie  must  be  initialized  (Table  1).  All  parameters  characterizing  the 
antenna,  the  subsurface  medium,  and  the  buried  objects  are  defined  in  this  file.  The  value  of  these  parameters  are  eas¬ 
ily  changed  without  altering  any  of  the  simulation  code.  The  definitions  of  each  parameter  follow. 
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Table  1 


Configuration  File  Parameters 


Sit*  Characterization  Parameters  (See  Figure  14) 

Parameter 

Meaning 

xbounds 

Bounds  in  the  x  direction  -  low  limit,  high  limit,  distance  step 

ybounds 

Bounds  in  the  y  direction  -  low  limit,  high  limit,  distance  step 

zbounds 

Bounds  in  the  z  direction  -  low  limit,  high  limit,  distance  step 

Antenna  Characterization  Parameters 

Parameter 

Meaning 

aheight 

antenna  height  (mm) 

times  tep 

time  step  (picoseconds) 

distbwant 

distance  between  antennas  (mm) 

anglestep 

angle  step  (radians) 

beamx 

beam  angle  -  x  cross  section  (radians) 

beamy 

beam  angle  •  y  cross  section  (radians) 

Reflection  Model  Characterization  Parameter 

Parameter 

Meaning 

speciobeangle 

specular  lobe  angle  (radians) 

Subsurface  Object  Characterization  Parameters 

Parameter 

Subparameters 

sphere 

center,  radius 

cone 

center  of  end  1,  radius  of  endl, 
center  of  end  2,  radius  of  end  2, 
length 

polygon 

number  of  points  (N) 
point  1  through  point  N 

patch 

number  of  points  (N) 
point  1 ,  normal  to  point  1 

point  N,  normal  to  point  N 

Subsurface  Characterization  Parameter 

Parameter 

Meaning 

baseperm 

base  permittivity 
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6.0  Extensions  for  Heterogeneous  medium 

When  a  volume  of  GPR  data  contains  subvolumes  with  different  electrical  characteristics  (heterogeneous  media), 
accuracy  of  the  processing  operations  can  be  greatly  increased  by  utilizing  a  more  complete  model  of  the  subsurface. 
The  electrical  characteristics  of  the  subsurface  are  used  to  calculate  the  velocity  of  the  radar  wave  in  the  medium. 
These  characteristics  are  used  in  the  heterogeneous  migration  algorithm  (Fenner  198S,  Ulriksen  1982). 

H  =  magnetic  susceptibility 
a  =  conductivity 

P  =  imaginary  part  of  the  wave  number 


A  subsurface  interface  is  a  division  between  two  subsurface  volumes  with  different  electrical  characteristics.  When  a 
radar  wave  is  incident  on  a  subsurface  interface,  the  amount  of  radar  energy  reflected  and  refracted  at  this  interface  is 
calculated  using  the  impedance  of  both  mediums  as  shown  below. 


Z,  =  Impedance  for  Layer  1  Z  = 

Zj  =  Impedance  for  Layer  2 

(Z,-Z,)  , 

Reflected  Energy  =  Incident  Energy 

Refracted  energy  =  Total  Energy  -  Reflected  Energy 


6.1  Forward  Modelling 

For  the  case  of  homogeneous  media,  when  a  ray  intersects  an  object  the  reflecting  ray  and  the  antenna  center  ray  are 
calculated  (Figure  16).  If  the  reflecting  ray  intersects  the  send  antenna’s  aperture,  then  the  energy  of  the  send  antenna 
is  placed  in  the  data  array  as  previously  described.  The  changes  for  heterogeneous  medium  area  conceptually  straight 
forward  (Figure  17).  Ray  intersections  are  done  recursively.  At  each  interface  two  rays  are  calculated:  the  reflected 
and  the  refracted  rays  (using  the  mode'  ^escribed  in  section  2).  There  are  two  end  conditions  for  the  recursion:  a  ray 
intersect  the  send  antenna’s  aperture  and  the  angle  p  is  within  tiie  user  specified  limits  or  a  ray  intersects  the  z-planc 
of  the  send  antenna  without  intersecting  the  aperture.  A  model  of  the  elecuical  parameters  for  each  subsurface  layer 
will  also  be  needed  since  the  index  of  refraction  (used  to  calculate  the  angles  of  reflection  and  refraction)  is  a  function 
of  the  medium’s  permittivity  and  the  conductivity. 
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FIGURE  16 


Forward  Modelling  In  Heterogeneous  Medium 
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6.2  Migration 

The  template  for  migration  in  homogeneous  media  consists  of  a  series  of  arcs;  every  point  on  each  arc  is  equidistant 
from  the  ground  surface’s  data  element  (Figure  17).  This  is  true  when  the  velocity  of  the  radar  wav.s  is  constant 
throughout  the  medium,  but  it  is  not  a  valid  assumption  for  the  case  of  heterogeneous  medium.  The  template  lakes  on 
some  arbitrary,  perhaps  irregular,  geometry  such  that  eve'y  point  on  each  arc  is  the  same  time  delay  from  the  corre¬ 
sponding  surface  data  element.  Ray  tracing  will  be  used  to  determine  the  coordinates  for  this  set  of  pomts.  In  the 
homogeneous  case,  a  template  was  calculated  in  a  preprocessing  step,  therefore  saving  computational  lime.  Because 
an  heterogeneous  algorithm  needs  to  work  for  any  distribution  of  electrical  parameters  in  the  subsurface,  a  prepro¬ 
cessing  step  is  not  possible.  This  means  the  template  must  be  calculated  for  each  data  element  corresponding  to  the 
surf^e  of  the  ground. 
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RGURE  17  Migration  Template  for  Heterogeneous  Medium 
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7.0  Discussion  and  Recommendations 

7.1  Discussion 

Since  the  current  focus  of  this  research  1  is  subsurface  mapping,  the  ability  to  accurately  locate  the  boundaries  of  bur¬ 
ied  waste  is  crucial.  The  test  site  used  for  this  evaluation  is  a  sandbox  with  the  dimensions  -  2  by  2  by  1.2  meters.  A 
cylinder  of  length  700  millimeters  and  radius  175  millimeters  was  buried  with  one  of  its  ends  coinciding  with  one  end 
of  the  .sandbox  (Figure  1).  The  left  and  right  edges  of  the  cylinder  were  predicted  (using  a  GPR  image)  with  less  than 
5%  error.  A  vertical  slice  of  unprocessed  data  taken  from  the  midpoint  of  the  cylinder  is  shown  in  Figure  2.  Data  from 
the  same  location,  after  deconvolution  and  three  dimensional  migration,  is  shown  in  Figure  3.  Unprocessed  data  for  a 
forward  modelled  sphere  with  the  same  radius  is  shown  in  Figure  4.  After  three  dimensional  migration,  the  same  data 
slice  is  shown  in  Figure  5.  The  coordinates  of  the  center  of  the  sphere  are  (1000,  600). 

A  significant  contribution  of  this  research  is  a  methodology  for  implementing  three  dimensional  migration  and  for¬ 
ward  modelling  of  GPR  data.  This  contribution  is  divided  into  three  areas;  accurate  modelling  of  the  antenna  beam’s 
attenuation  pattern,  accurate  modelling  of  the  subsurface,  and  computationally  fast  ray  tracing. 

Digitized  radar  reflections  in  actual  GPR  data  consist  of  positive  and  negative  values.  The  migration  operation  causes 
positive  reinforcement  in  the  area  where  reflections  due  to  subsurface  objects  or  interfaces  occiirrcd.  The  operation 
causes  negative  reinforcement  in  the  areas  where  there  arc  no  reflections  due  to  subsurface  objects  or  interfaces.  Of 
course,  in  areas  where  there  are  only  weak  reflections  the  reinforcement  caused  by  the  migration  operation  is  not 
complete.  Al!  the  data  generated  by  the  forward  modelling  algorithm  will  have  positive  values.  When  the  migration 
operation  is  used  on  data  produced  by  the  forward  modelling  algorithm,  the  positive  reinforcement  works  but  the  neg¬ 
ative  reinforcement  does  not. 

7.2  Recommendations 

Implementation  of  three  dimensional  heterogeneous  migration  and  forward  modelling  is  the  next  important  step.  iLs 
algorithm  is  outlined  in  the  previous  section.  Because  the  velocity  of  the  radar  waves  vary  significantly  in  different 
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medium,  the  accuracy  of  the  located  subsurface  objects  should  significantly  improve.  The  path  of  the  traced  rays  will 
change  significantly  at  many  interfaces.  Since  the  accurate  location  of  subsurface  objects  is  of  utmost  importance,  a 
quantitative  method  is  needed  to  measure  this  accuracy.  A  computerized  method  is  needed  to  compare  the  actual 
location  and  geometry  of  the  object  (apriori  knowledge  assumed)  with  the  location  and  geometry  obtained  from  the 
GPR  images.  This  could  be  done  in  three  steps.  Store  the  actual  location  and  geometry  of  each  object  in  a  file,  display 
it  in  an  overlay  on  the  GPR  image  (computer  screen),  and  use  visual  inspection  to  determine  the  image’s  accuracy.  A 
better  method  is  to  eliminate  visual  inspection  by  using  an  edge  detection  algorithm  on  the  GPR  image  before  com¬ 
paring  it  with  the  computer  model. 
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